Update ImPACT code to work again #2305
Conversation
…, introduced dummy interpolators to use when testing so templates are not needed.
…nto ImPACT_tests
…nto ImPACT_tests
| * template_network_interpolator - Now works with templates with different zenith and azimuth angles | ||
| * unstructured_interpolator - Significant speed improvements | ||
| * pixel_likelihood - Constants added back to neg_log_likelihood_approx, these are quite important to obtaining a well normalised goodness of fit. | ||
| * hillas_intersection - Fixed bug in core position being incorrectly calculated, fixed tests too |
There was a problem hiding this comment.
Could you add a sentence about atmosphere_profile now being an argument to Reconstructor()?
| default_value=False, help="Use time gradient in ImPACT reconstruction" | ||
| ).tag(config=True) | ||
|
|
||
| root_dir = traits.Unicode( |
There was a problem hiding this comment.
Yes. Maybe even removing the hard coded names for the templates is worth a look
There was a problem hiding this comment.
That could be a TelescopeParameter(Path)
| root_dir = traits.Unicode( | ||
| default_value=".", help="Directory containing ImPACT tables" | ||
| ).tag(config=True) | ||
|
|
||
| # For likelihood calculation we need the with of the | ||
| # pedestal distribution for each pixel | ||
| # currently this is not available from the calibration, | ||
| # so for now lets hard code it in a dict | ||
| ped_table = { |
There was a problem hiding this comment.
This should at least be a FloatTelescopeParameter
There was a problem hiding this comment.
And further, it woudl be best to first try to get the peds from the event (i.e. from event.mon.tel[1].pedestal.charge_std), and only fall back to the ped_table if not.
These values are probably not even correct for Prod6.
| template_scale=1.0, | ||
| xmax_offset=0, | ||
| use_time_gradient=False, | ||
| self, subarray, atmosphere_profile, dummy_reconstructor=False, **kwargs |
There was a problem hiding this comment.
I think it's relevant for a couple of tests
There was a problem hiding this comment.
It seems to be a bit misnamed, as it replaces the TemplateNetworkInterpolator with a dummy one with a dummy template, so maybe call it "no_template_mode" or "test_mode" or something like that.
| # self.priors = prior | ||
|
|
||
| # String templates for loading ImPACT templates | ||
| self.amplitude_template = Template("${base}/${camera}.template.gz") |
There was a problem hiding this comment.
Is templating only by camera name what we want? Shouldn't we at least use the full str(telescope_description)?
E.g. FlashCam on MST vs. FlashCam on HESS II
| telescope_pointings = self._get_telescope_pointings(event) | ||
|
|
||
| # Finally get the telescope images and and the selection masks | ||
| mask_dict, image_dict, time_dict = {}, {}, {} |
There was a problem hiding this comment.
In the HillasRecontructor we gained large speed improvements by not using dicsts tel_id -> value but dense arrays, maybe could also help here?
| if mean_height > 100000 or np.isnan(mean_height): | ||
| mean_height = 100000 | ||
| if mean_height_asl > 100000 or np.isnan(mean_height_asl): | ||
| mean_height_asl = 100000 |
There was a problem hiding this comment.
why replace nan or nonsensical values with a valid looking value? That might be misleading. Also a magic number.
| pass | ||
|
|
||
|
|
||
| def guess_shower_depth(energy): |
There was a problem hiding this comment.
Astropy units? Or is this runtime critical?
This should maybe better be an internal function.
Is there a reference for this hard coded values?
| # If our energy estimate falls outside of the range of our templates set it to | ||
| # the edge | ||
| if lower_en_limit < 0.02: | ||
| lower_en_limit = 0.02 |
There was a problem hiding this comment.
Shouldn't this be somehow stored in the template instead of being hardcoded?
There was a problem hiding this comment.
E.g. if I produce templates for SSTs, lower energy of 20 GeV really doesn't make sense...
| x_max_exp = guess_shower_depth(energy) | ||
| diff = xmax - x_max_exp | ||
| return -2 * np.log(norm.pdf(diff / width)) | ||
| # These are settings for the iminuit minimizer |
There was a problem hiding this comment.
I assume these are all fixed settings, i.e. things you don't ever expect the user to adjust?
If that is not the case, then you should add them to __all__ and change the comments to use the form #: Comment so they become documented.
| root_dir = traits.Unicode( | ||
| default_value=".", help="Directory containing ImPACT tables" | ||
| ).tag(config=True) | ||
|
|
||
| # For likelihood calculation we need the with of the | ||
| # pedestal distribution for each pixel | ||
| # currently this is not available from the calibration, | ||
| # so for now lets hard code it in a dict | ||
| ped_table = { |
There was a problem hiding this comment.
And further, it woudl be best to first try to get the peds from the event (i.e. from event.mon.tel[1].pedestal.charge_std), and only fall back to the ped_table if not.
These values are probably not even correct for Prod6.
| template_scale=1.0, | ||
| xmax_offset=0, | ||
| use_time_gradient=False, | ||
| self, subarray, atmosphere_profile, dummy_reconstructor=False, **kwargs |
There was a problem hiding this comment.
It seems to be a bit misnamed, as it replaces the TemplateNetworkInterpolator with a dummy one with a dummy template, so maybe call it "no_template_mode" or "test_mode" or something like that.
|
Hi all, I am going to merge this now, thanks to everyone who put in the hard work to make this function again! @gschwefer @Tobychev @ParsonsRD There are a lot of comments for improvements here, but these are very hard to review, so I will merge here and let's continue in easier to review follow-up PRs. |
As requested by @kosack I open a PR with my version of @ParsonsRD ctapipe ImPACT implementation, where I've introduced changes to make the code able to merge with the current HEAD.
This branch also has a modification to some reco tests, and in particular
test_reconstruction_changes.pyexpects that you have a local copy ofLSTCam.template.gzand NectarCam.template.gz that I've generated from a larger template file supplied by @ParsonsRD, both should now be available on the test server.The events used for the test are very low energy events, and the NectarCam template is actually just the LST one again, so performance on these events is not stellar, but they are processes without crashes.